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The developing structure in systems of compacting ductile grains were studied experimentally in 
two and three dimensions. In both dimensions, the peaks of the radial distribution function were 
reduced, broadened, and shifted compared with those observed in hard disk- and sphere systems. 
The geometrical three-grain configurations contributing to the second peak in the radial distribu- 
tion function showed few but interesting differences between the initial and final stages of the two 
dimensional compaction. The evolution of the average coordination number as function of packing 
fraction is compared with other experimental and numerical results from the literature. We con- 
clude that compaction history is important for the evolution of the structure of compacting granular 
systems. 
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I. INTRODUCTION 

The structure of mono-disperse granular media is 
known to be very sensitive to the shape of the grains 
0, grain interactions such as cohesion 0, and as- 
sembling procedure |^. Packing fractions as low as 
c = 0.125 have been found experimentally in 3D-systems 
with high grain-grain attraction (5|, whereas dense sys- 
tems of spherical grains can be packed to c > 0.65 by 
a carefully selected tapping procedure. 

Historically, the structure of dense granular media was 
studied during the 1960s and 70s as a model of fluids 
and amorphous materials 0, 0- Since then, the com- 
plex properties of granular structures have been studied 
in dense elastic pack ings with much focus on compaction 
dynamics 0, 0, The densification of granular 

packings has a variety of applications, and in particular, 
it takes place in Nature during the slow compaction of 
sediments 12] and the fast event of a landslide. Gran- 
ular compaction is also commonly studied in relation to 
pharmaceutical powders and metal industry |l3|. In the 
compaction of sediments and powders, only the initial 
stages can be modeled by the compaction of hard elas- 
tic grain ensembles. When the geometrical structure is 
jammed |14| (i.e., no grain can be geometrically trans- 
lated while all others remain fixed), at packing fractions 
c ^ 0.64, further compaction can only occur by grain de- 
formation (l^j . For the compaction of sediments and the 
hot isostatic pressing of metal powders, grains deform 
plastically, thus deformed regions are relaxed during the 
compaction. 

The development of a dense structure in a ductile gran- 
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ular ensemble will depend on grain properties such as 
friction and cohesion, as well as the deformation pro- 
cedure. In such ductile packings, the densest possible 
structure will depend on how pore space can be reduced 
at high packing fractions. At such packing fractions, un- 
less the granular ensemble is compacted in vacuum, gas 
or fluid occupying the pores must be transported out or 
absorbed into the matrix for further compaction to oc- 
cur. An example of such a densiflcation experiment is 
the experiment by Stephen Hales in the early 1700s [T5| . 
in which he filled a glass with dried peas and water. The 
peas absorbed water and eventually filled the container 
as an ensemble of polyhedra. Very little is known about 
the structure of plastically deforming grain ensembles, 
although they are extremely important for a range of in- 
dustrially and naturally occurring compaction processes. 



Here, we present a series of experiments on ensembles 
of compacting ductile grains in two and three dimensions. 
Positions and coordination numbers have been measured, 
and structural measures such as the radial distribution 
function and distribution of coordination numbers were 
obtained. In particular, we were able to study the lo- 
cal geometry that dominates the peaks in the radial dis- 
tribution function, and the evolution of this geometry 
with packing fraction in two dimensions. The results are 
compared to similar findings from the literature of both 
hard and ductile experimental grain ensembles, and to 
a numerical compaction model (Rampage) , developed to 
simulate the initial compaction of sediments . 



The following two sections ljIIIIII|l presents the two and 
three dimensional systems. The discussion then follows 
in section Hvl 



II. THE TWO-DIMENSIONAL SYSTEM 

The two-dimensional (2D) system consisted of 1100 
oriented cylinders, which were compacted at constant 
rate. From image analysis, the positions, coordination 
numbers, and packing fractions were found at regular in- 
tervals during the compaction. 

A. Experiment 

The 2D-setup consisted of ductile cylindrical grains, 
which were stacked in a Hele-Shaw cell, and uniaxially 
compressed by a piston at a constant rate v = 1 mm/h. 
The grains were made of spaghetti, which had been boiled 
in water for 14.5 minutes in order to render them soft 
and deformable. The grains' resistance to deformation 
could be controlled by the boiling time. After boiling, 
the spaghettis were kept in a clear olive oil to avoid de- 
hydration, then cut into cylinders of length 10 mm, and 
carefully stacked in the container (Hele-Shaw cell). Two 
sizes of grains were used, but in different sections of the 
model. For simplicity, only the lower section of nearly 
600 mono-disperse grains was used in this analysis. These 
grains had (initial) diameters of d = (2.8 ± 0.1) mm. 

The width and height of the initial packing was w — 83 
mm and h = 92 mm, corresponding to 29 and 33 grain di- 
ameters, respectively. The initial height of the 600 lower 
grains was 55 mm, or 20 grain diameters. During the 
stacking, olive oil was constantly added to the container 
so that the grains were immersed in oil at all times. The 
oil, apart from preventing dehydration of the grains, also 
lubricated the walls so that friction did not restrain the 
compaction. The depth of the Hele-Shaw cell was 11-12 
mm, to ensure a channel at the back of the stacked grains 
for the oil to escape through. The oil thus did not affect 
the compaction mechanically. Figure ^a) illustrates the 
setup of the 2D experiment, and Fig. [21contains a closeup 
of a region in the first image. Pictures were taken every 
20 minutes, and a total of 47 pictures were taken as the 
system compacted at a speed of w = 1 mm/hour from an 
initial packing fraction of c = 0.89 to the final c = 0.99. 
The camera was an AstroCam Capella, LSR Life Science 
Resources, UK, with 3000x2000 pixels, 14 bit resolution. 
The spatial resolution was 23 pixels per mm. 

B. Analysis (2D) 

The geometrical center position was obtained for all 
grains by image analysis. In the first image, i.e., the 
image taken immediately after compaction started, the 
center of mass was found by the following procedure: 
considering circles around a particular pixel, the radius 
of the circle that contained 5 pixels of intensity below a 
certain threshold was taken to be the shortest distance 
to the pore space from this particular pixel. The pixel 
within a grain with the largest distance to the pore space 
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FIG. 1: Schematic illustrations of the experimental setups, h 
is the filling height, and d is the grain diameter in both 2D 
and 3D. (a) The quasi two-dimensional arrangement, w is the 
width of the container, and the piston is driven at constant 
velocity v. The grains were cylinders made from spaghetti, 
and stacked such that their length was in the depth of the 
container. When viewed from the front of the container only 
the circular crossection of the grains were seen, of diameter 
d. (b) The three-dimensional cylindrical ensemble, consisting 
of spherical grains of diameter d made of Playdoh. The inner 
container diameter is D. 

was taken to be the center position of the first image, see 
Fig. [3 for an illustration. None of the disks were much 
deformed at this stage, so looking for pore space in cir- 
cular sections accurately determined the center position 
to a resolution of a fraction of a pixel for most grains. 
However, due to reflected light from the Hele-Shaw cell, 
some positions were only determined within a few pix- 
els resolution. The center of mass positions found in the 
first image were traced by pattern recognition (normal- 
ized cross correlation function in the following 
images, to a resolution of 1 pixel. 

The grains were essentially incompressible. During 
compaction, the average crossectional area A of each 
grain decreased 6%, thus the grains elongated in the di- 
rection perpendicular to the image plane. The develop- 
ment of the crossectional area was found manually for 
ten specific grains in the first and last image. Of these 
ten grains, two representative grains were selected, and 
their areas were found (also manually) in every second 
image throughout the experiment. A best fit was made 
to the developing average area, as A,n — a(l — {m/b)^), 
where m is image number, a — 3.6 10^ square pixels, and 
b — 195. The effective grain diameter d{m) in each image 
was then found as d{m) = a/ 4Am/7r. Figure 13 shows the 
grain diameter, calculated from the average area of each 
image, and the effective diameter ^(to) (line) as functions 
of image number m. By using the effective diameter in 
all the analysis, the decreasing crossectional area does 
not affect the geometry and structural evolution during 
the compaction. Thus, the packing fractions and radial 
distribution functions can be compared between images 
independent of the actual grain diameters. 
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FIG. 2: Closeup of a grain in the first image. The procedure 
of finding the grain center in the first image is illustrated: the 
two circles represent the largest distance from their center 
pixels at which a maximum of 5 "dark" pixels are found (see 
text for details). 



1. Packing fractions 

The packing fraction c of the 2D system was obtained 
by Voronoi analysis of the position data from each image. 
Only disks at a distance of more than one grain diam- 
eter from any boundary were used in the calculation to 
avoid unbounded 18] Voronoi cells. The area A^i of the 
Voronoi cell for each disk center i was found, and as each 
such cell contains one disk, the local packing fraction Ci 
was given by 

Ci^iTcP/AA^„ (1) 



where d is the average diameter of a grain in the image 
in question. The Voronoi tessellation of a region around 
a certain grain is shown for the first and last (dashed) 
image in Fig. E^a), with grain centers marked as bullets 
(first image) and circles (last image). The motion rel- 
ative to the central grain is indicated by lines between 
the grain center positions in the first and last image. 
Figure Elb) shows the cumulative distribution P{r) of 
distances between touching grains for the first (curve A) 
and last (curve B) images. The overall packing fraction 
based on the Voronoi cell division was c = (q). The 2D 
compaction spanned packing fractions from c — 0.89 to 
c = 0.99. 
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FIG. 3: The effective grain diameter, d{m), (line) as func- 
tion of image number, m. (•) Grain diameters, d = \/ AA/tt, 
calculated from the average of the measured area. A, of two 
grains. 




FIG. 4: (a) An example of the Voronoi structure in a region of 
the first and last (dashed) images, used in obtaining the local 
packing fractions d. The tessellation of the last image was 
translated so that the center of one cell lie undistorted. The 
compaction direction is downward. The cumulative distribu- 
tions of touching neighbor distances are shown in (b), for the 
first (A) and last (B, dashed) images. The average value of 
distances between touching neighbors is indicated by vertical 
lines for the first and last (dashed) image. 



2. Radial distribution function 

The radial distribution function (RDF) g{r) is defined 
by g{r) = p{r)/{p), where (p) = N/V is the average 
number density of grain centers in a container of area V, 
and p{r) is the average number density as a function of 
distance r from a grain center. The expression for the 
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number density as measured from grain number i is 



Vn{r) 



(2) 
(3) 



Here, Vn{r) — 2'Kr is the circumference of a circle of 
radius r, and Ni{r) is the number of grain centers at 
distance r from grain i. r.ij is the distance between grains 
i and j, and 5 is the Kronecker delta, 5=1 for rij = r, 
and otherwise zero. The average Pi{r) over N grains gives 



p{r) 



1 ^ 



thus 



1 



N{p)Vn{r) 



(4) 



(5) 



where i is summed over all grain centers of the sample. In 
a finite system the expression for g{r) should be replaced 
by 



l{r) = 



N{p)Vn{r, dr) 



^ ^N^{r (^[r- dr/2, r + dr/2]), 

i 

(6) 

where dr is the width of the shell, its size depending 
on the number of grains in the system and/or the uncer- 
tainty of the position measurements. Perturbations to dr 
should not affect g{r) when the proper value of dr is cho- 
sen, as the dependence of Ni and Vn on dr should cancel 
in the expression of g{r). The expression for the circular 
shell becomes Vn{r, dr) = 7r[(r + ^ {f ^ = 

2'Kr dr. 

Noting that the sum in Eq. ^ over Ni equals twice the 
number of distances n(r, dr) of lengths r G [r — dr/2, r + 
dr/2] in the ensemble, X); '^^) — 2n{r,dr), the ex- 
pression for the RDF of a finite packing becomes 



g{r) = 2n{r,dr)F{r)/{N{p)^Vn{r,dr)) 



(7) 



where F(r) is a finite size correction for the boundaries, 
as discussed below. 

When Eq. {Tj) is used without the normalization F{r), 
the boundaries of the ensembles introduce finite size ef- 
fects to the RDF of a small system. These finite size 
effects result not only because the structure along the 
boundaries differ from the interior structure, but because 
the normalization over circular shells of areas 2'Krdr in- 
cludes regions outside of the ensemble. To avoid the lat- 
ter source of errors, the normalization function F{r) is 
introduced based on the specific rectangular geometry 
of the ensemble. F{r) is found as follows: In an infi- 
nite system, a grain in a position i? = {x,y) would be 
surrounded by a circular shell 2'Kr dr independent of its 
position. In the finite system of a container, whether all 



of the circular shell lies within the container depends on 
the radius r and the grain's position R. The fraction of 
circular shells of radius r that partly lie outside of the 
container increases with r and is a measure of the error 
done by disregarding the boundaries. Let the area of the 
container be V . The integral over V of all circular shell 
circumferences 2Kr with center positions inside the area 
of the container is 2'KrV . Only a part of this integral rep- 
resents circumferences that lie inside the container, thus 
in normalizing the RDF of a finite system one should use 
this fraction instead of the 2KrV . A normalization func- 
tion Fir) can be defined as the ratio of the integrals of 
circumferences in the infinite case to the finite case: 

Sy2'Kr d'^R 



F{r) 



Jy A{r, R) d^R 



(8) 



A{r, R) is the fraction of the circular shell of radius r cen- 
tered at R that lies inside the container. For a rectangu- 
lar container of width tu and height h the normalization 
function is: 

F{r) = Kwh[K{w - 2r){h - 2r)+{K - l){w + h - Ar) 

+ r(37r/2-2)]"\ 

(9) 

This correction procedure was previously used by Mason 
[T^ for a different geometry (3D spherical ensemble). 



3. Coordination numbers 

The coordination numbers k (number of contacts per 
grain) were found from further image analysis; a Delau- 
nay triangulation based on the position data was used to 
obtain the set of nearest neighbors of each grain. Each 
Delaunay neighbor connection was then investigated by 
a thresholding procedure to establish whether it was a 
touching neighbor: The intensity values of the array of 
pixels that formed the shortest path between grain cen- 
ters i and j showed a distinct dip whenever pore space 
was present between the grains. Grains were considered 
not touching if the dip was below 2/3 of the intensity 
value of the grain face intensity. Figure illustrates this 
procedure. Grains i and j, as seen in Fig. |3Ja), are not 
touching as the minimum value of their center-to-center 
intensity plot is below I = 850-2/3 = 567, as seen in (b). 
The intensity values was in the range / S [100,1000], 
where / = 100 represented black. The typical difference 
in intensity values between open pore space and the inter- 
nal of a grain was A/ = 500. The coordination number 
distribution P(fc) was found with this procedure for all 
the images. 



C. Results (2D) 

The radial distribution functions for the first and last 
stages of the experiment are shown in Fig. |3 The struc- 
ture at c = 0.89 shows many of the characteristic peaks 
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FIG. 5; The procedure of finding touching neighbors was 
based on the intensity values of the direct line between grain 
centers, (a) shows a small region of the initial compaction 
state, where i, j, and k are grain centers. The lines i-j and 
i-k are the direct lines between the centers, of which the in- 
tensities are plotted in (b) as curve A and B, respectively. 
The dip shown in curve B was typical when grains did not 
touch. 



of a triangular packing of disks, except that the peaks are 
broadened and shifted to shorter distances than those of a 
crystalline packing. Crystalline correlations occurring in 
triangular lattices of hard disks at r/d = (1, a/3, 2, V?, 3) 
are shifted to r/d = (0.958, 1.695, 1.915, 2.580, 2.875). In 
the RDF of the more compacted system of c = 0.99 the 
peak near r/d = 3 has vanished, and the double peak 
near ~ 1.7 has become smoother and broader. Also 
the first peak {r/d = 1) is seen to broaden considerably 
as many neighbor distances at this stage are less than 
one effective grain diameter. The inset of Fig. shows 
the finite size correction function F{r) for the first image, 
which is practically identical to the correction in the last 
image. 

The geometrical structures associated with the peaks 
at r/d = 1.695 and r/d = 1.915 in the RDF are illus- 
trated in Fig. by a conditional three-point correlation 
function Figure [T^a) was generated by identifying 

all grain centers at distances r/d E [1.695 ±0.05] to each 
other in the packing, aligning each pair of such centers 
along the horizontal axis, and plotting the surrounding 
grain centers in the corresponding positions. The two 
fixed grain positions are shown on the horizontal axis at 
the center of the figure, while the surrounding intensity 
structure visualizes the probability of having a third grain 
in any position relative to the two fixed ones. The inten- 
sity was normalized by the number of pairs contributing 
to the plot. Figure [7fb) is generated in the same manner 
as (a), but for grains at a distance of r/d G [1.915±0.05] 
to each other. The upper half of (a) and (b) shows the 
probability distribution for the initial stage of the com- 
paction (c = 0.89), while the lower part shows the distri- 
bution at c = 0.99. Black represents a high probability 
of the third grain, and white represents zero probabil- 




FIG. 6: The radial distribution function, g{r), as a func- 
tion of particle separation r in the ductile 2D experiment, 
for packing fractions c = 0.89 (o) and c — 0.99 (•, bold 
line) . The curves connecting the data points are produced by 
splines. The normalization function for the 2D setup, F{r), 
Eq. is shown in the upper right corner (curve A) for 

the distances r/d £ [0,5]. The typical geometrical configura- 
tions that contribute to the peaks at r/d £ [1.695 ± 0.05] and 
r/d e [1.915 ± 0.05] are illustrated. 
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FIG. 7: Visual illustration of the probability of a third grain 
position for two grains at a given distance (a) r =1.645-1.745 
d, and (b) r =1.865-1.965 d. The upper part of the plots 
show the probability distribution at c = 0.89, while the lower 
give the distributions at c = 0.99. 



ity. The black circles mark the distance of one diameter 
from the two fixed grain positions. The intensity pattern 
suggests a regular structure with obvious symmetries in 
the packing. At first glance, the structure might seem 
hexagonal, but this is due to the alignment of the two 
fixed centers; The underlying structure is triangular. Al- 
though the intensity plot of Fig. [71[a) does not prove the 
presence of the four-point configuration illustrated in Fig. 
closer inspection of the configurations contributing to 
this peak (at r/d = 1.695) in the RDF shows that this 
configuration is dominant for r/d — 1.695. Similarly, the 
configuration of three grains in a row is seen to contribute 
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FIG. 8: The distribution P{k) of coordination numbers k as 
function of packing fraction c in two dimensions. Curves A 
and B are projections of P{k) at c = 0.89 and 0.99, respec- 
tively. Curves C, for P{5), and D, for P{6), are developments 
of P{5) and P(6) with increasing c. 
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FIG. 9: Mean coordination numbers (A;) as a function of pack- 
ing fraction c for the two dimensional ductile experiment. The 
circle represents the 2D space filling structure. 



most to the peak at r/d = 1.915 of the RDF, although 
the center grain is often slightly misplaced along the axis. 
The lower part of Fig. CJa) contains a small probability 
of a third grain in between the two fixed grains (barely 
visible in the figure), which is not seen in the upper part 
(c — 0.89). This very compressed alignment of three 
grains in a row emphasize the role of ductility of the 
grains, as this configuration is impossible in packings of 
hard disks. 

The cumulative distribution P{r) of distances between 
touching grains is shown in Fig. Elb). The average dis- 
tance between touching grains decreases from the first 
to the last image, while the width of the distribution in- 
creases. The touching grains at distances larger than a 
grain diameter are grains aligned perpendicular to the 
compaction direction, thus the ellipticity of grains in- 
creases, especially in the last compaction stage. 

The coordination number distribution is shown in Fig. 
IHl as a three-dimensional plot of the distribution P{k) 
as function of coordination numbers k and packing frac- 
tion c. The distribution is broad at c = 0.89, and 
then narrows and shifts to higher coordination numbers 
with increasing packing fraction. At packing fractions 
c = 0.906-0.918 a few cases of grains with only two con- 
tacts are observed, but only one or two such grains are 
present at each packing fraction. Curves A and B in Fig. 
IHlare projections of the first and last distributions of the 
experiment, respectively. Curves C and D are projections 
of the fractions of five and six coordination numbers with 
increasing packing fraction. 

Figure^shows the average coordination number {k) as 
function of packing fraction in the two dimensional exper- 
iment. The circle at coordinates (1,6) represents the end- 
point for 2D systems, as a space filling structure (Voronoi 
cells) in 2D must have six neighbors on average [2ll |. 



Also shown are the data for two crystalline structures, 
the square lattice and the triangular lattice. We observe 
that the evolution of the average coordination number in 
the ductile packing closely follows a straight line between 
the square lattice and the space filling Voronoi structure. 



III. THE THREE-DIMENSIONAL SYSTEMS 

The three-dimensional (3D) systems analyzed here are 
six experimental systems of compacting ductile grain en- 
sembles. These are compared to Finney's experiment 
7] on a random close packed ensemble of hard spheres 
(c = 0.6366), from which the positions were kindly made 
available to us. This ensemble contained nearly 8000 
grains in a square box, and the positions of the grains 
were measured to a precision of 0.2% of the grain diam- 
eter. A numerical model of a compacting grain ensemble 
is also presented, giving possible trajectories for com- 
paction between packing fractions c = (0.50-0.64). 



A. Experiments 

The 3D experimental setup consisted of N ductile 
grains poured into a cylindrical Plexiglas container, ap- 
proximately ten at a time, to a height h. The set of 
grains were then left to compact under gravity for a time 
t, or were mechanically compacted to a specific packing 
fraction. A lid covered the top of the cylinder during 
the experiment to avoid dehydration, as the grains con- 
tain water and other volatile chemicals. Figure ^^b) il- 
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lustrates the setup of the 3D experiments. The grains 
were prepared from Play-Doh (Hasbro Int. Inc., UK) to 
spheres of diameters d — (8.8 ±0.2) mm. The material is 
viscous for shear rates in the range (10~'^-10~^)/s, with 
a viscosity of 3 • 10^ Pas, as measured in rotational tests 
(Paar-Physica MCR 300, rheometer). Except from the 
first few minutes and the very end of an experiment, the 
strain rate was within this regime. Strain rates down 
to 10~''/s were common after a week's compaction, but 
shear tests could not be performed at such low rates, and 
the material properties for these rates are not known. 

Each experimental system was disassembled grain by 
grain while measuring positions and/or coordination 
numbers of the ensemble. The position coordinates 
{x,y,z) of the grains' top point were established with 
a mechanical arm (MicroScribe) to an estimated resolu- 
tion of 0.5 mm (6% of grain diameter), limited by the 
difhculty of establishing the center of the grain's top sur- 
face. The number of touching neighbors (coordination 
number) was also counted for each grain during the disas- 
semblage. To ease this procedure, grains of five different 
colors were used, as grains of different colors left marks 
on each other in contacting areas. Before a grain was 
removed, the visibly contacting neighbors were counted. 
The grain was then removed, and any formerly hidden 
contacts became visible due to grain deformation and/or 
dis-coloring. As a grain was removed, all its contact- 
ing grains that remained in the ensemble were carefully 
marked with a felt pen. Thus, the coordination number 
of a grain consisted of the number of marks on its upper 
surface, plus the number of visible and hidden contacts 
established during the removal of the grain itself. The 
number of contacts that each grain had to the cylinder 
wall and/or floor was also recorded. 

Six different three dimensional ensembles were studied: 
Two of them with 2200 grains in a cylinder of diameter 
D ~ 130 mm, the rest with 150 grains in cylinders of 
diameter D = 32 mm. The four small systems were pre- 
pared simultaneously, and left to compact for different 
amounts of time before each was disassembled and the 
positions and coordination numbers of each system mea- 
sured. The resulting packing fractions were c =0.64, 0.64, 
0.66 and 0.68 ±0.01. The small system size was desirable 
in order to avoid considerable compaction during the dis- 
assemblage for the short compaction times. One of the 
big systems was left to compact for 145 hours before it 
was disassembled, measuring the positions and coordina- 
tion numbers of the grains. The packing fraction went 
from 0.57 (based on filling height) to 0.701 during the 
compaction in this ensemble. Due to a calibration error 
with the mechanical arm, only the positions of 787 grains 
in the mid to lower part of the cylinder were successfully 
measured, while the coordination number was established 
for 1169 grains. The second large system was compacted 
to a packing fraction of c = 0.75 by gently pushing a 
piston from the top of the cylinder. Compacting the sys- 
tem to such high packing fractions would otherwise have 
required months of gravitational compaction. Only the 



coordination numbers were measured in this ensemble, 
for 839 central grains. 



B. The Finney ensemble 

Finney carefully measured the positions of nearly 
8000 mono-disperse spheres in a box of square crossec- 
tion. The spheres were steel ball-bearings of diameters 
1/4 inch, which positions were measured to 0.2% of the 
grain diameter. The central 817 grains were used in this 
analysis, as these grains had a minimum distance of five 
grain diameters to any boundary, thus no corrections to 
the RDF was necessary to r/d = 5. Based on this se- 
lection, the radial distribution fimction and coordination 
numbers were found. 



C. Numerical compaction model 

A numerical model (Rampage, 16] ) was used to sim- 
ulate dry granular media. It consisted of randomly plac- 
ing one thousand mono-disperse spherical grains of di- 
ameter d in a container of square cross-section. Peri- 
odic boundary conditions applied in the horizontal plane. 
Grains were randomly positioned until a packing fraction 
of c = 0.5 was obtained. No considerations about over- 
lapping particles were done during this initial filling pro- 
cedure. The overlaps were then reduced by repositioning 
grains according to a force balance based on the elastic 
energy of the overlap region, and gravity. An overlap of 
volume Vo resulted in an inter-grain force / of 



in which Y is Young's modulus of the material. When 
gravitationally unstable, a grain would drop or roll un- 
til a gravitationally stable condition was obtained. This 
compaction process went on until a predefined tolera- 
tion threshold of average overlap distance I was reached 
{l/d < 0.005), while the ensemble was gravitation- 
ally stable. At this point further compaction was ob- 
tained by tapping, an incremental (< O.OOl/i) vertical 
re-positioning downward of each particle (at height h), 
followed by the same procedure of reducing overlaps and 
reaching gravitational equilibrium. Packing fractions of 
stable numerical ensembles were in the range c — 0.55- 
0.64. This model was originally developed for the mod- 
eling of sediment compaction, and is described in detail 
in Ref. 1^. The model reproduces statistically Finney's 
ensemble at c = 0.636 with respect to the RDF, the co- 
ordination number distribution and the distribution of 
contact angles, see Ref. |T^. 
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D. Analysis (3D) 

1. Packing fractions 

In three dimensions, the packing fractions were com- 
puted from the position measurements for all the cases 
in which these were available, which were most of the 
ductile experiments, the Finney ensemble, and the nu- 
merically generated ensembles. Boxes of variable square 
cross-section P and height h equal to the height of the 
packing were centered in the granular ensemble for cal- 
culation of the packing fraction. For each of these boxes 
j, the packing fraction was Cj{l) — NjVg/Vbj, where Nj 
is the number of grains contained in the box, Vg = ttcP /6 
is the average volume of a grain, and Vbj — Ph is the 
volume of the box. Grains that partially intruded the 
box, i.e., when their center position was less than one 
grain diameter from the box boundary, also contributed 
to the number of grains in the box. The overall packing 
fraction c of an ensemble was obtained by averaging over 
n — 10 boxes in distances 0.5-1.5 grain diameters d from 
the boundary: 



1 



(n + 1) 



Cj{lo+jAl). 



(11) 



where Iq is the side of the box in a distance 1.5 grain 
diameters from any boundary, and Al — d/n. In the nu- 
merical model, the total volume of two touching grains 
was smaller than two grain volumes if the grains over- 
lapped. The overlap volume was not assumed to be 
transported to the pore volume, as this would effectively 
change the structure of the packing, thus the calcula- 
tion of the packing fraction must account for the excess 
overlap volumes in the packing. This was done by reduc- 
ing the sphere radii by the mean overlap distance in the 
packing, thus reducing the volume of the grains before 
calculating c. 



2. Radial distribution function 

The radial distribution function in three dimensions is 
given by Eqs. © and Q, but with Vn as the volume of 
a spherical shell of radius r and width dr in 3D: Vn — 
47r[(r -t- dr/2)3 - (r - dr/2)3]/3 = 47rr2 dr -nidrf jZ. 

Equation ® was used in the calculation of the RDF in 
the Finney ensemble, as an ensemble average was made 
only over the central 817 grain positions to avoid bound- 
ary effects. These grains were all in a minimum distance 
of five grain diameters from any boundary, thus the RDF 
is not affected by boundaries up to r/d = 5. In or- 
der to calculate the RDFs of the large ductile ensemble 
(c = 0.70), Eq. lO was used, as the statistics of the whole 
ensemble was needed for a sufficiently detailed RDF to 
be found. A finite size correction function F(r) was in- 
troduced for the ductile ensemble in the same manner as 



for the 2D system, Eq. ©. In 3D, A{r,R) is the frac- 
tion of the spherical shell area (radius r, center position 
at i? = {x,y), height z) which is inside the cylindrical 
container. The resulting correction function is 



Fir) 



AiT^r^R^h 



Jo Ui 



h)dz 



where 



Il = 2n^r{h - z)R^ 



-2 / ATir{h - z){R-ux) 
lo 

'2xR - u{x'^ - 1) 



(12) 



(13) 



X arccos 



2{R-u) 



(14) 



dx. 



and u = y/r"^ — . F{r) is presented in the inset of 
Fig. El ^ is here the length of ii, i.e., the radial distance 
of a grain center from the cylinder axis. The width of the 
spherical shells were chosen so that the standard error of 
distribution from the mean was less than two per cent of 
the RDF, for the second peak. This gave the following 
resolutions of the RDFs; dr jd = 0.1 (ductile ensemble) 
and dr/d = 0.02 (Finney). 



3. Coordination numbers 

Coordination numbers in three dimensions were found 
experimentally for the ductile ensembles, as described in 
section IlII Al As the grains in contact to the walls were 
known, the distributions of coordination numbers for the 
ductile ensembles are based only on the internal grains. 
In the second of the large ensembles (c = 0.75), only 
grains at a minimum distance of three grain diameters 
from any boundary were included. For the Rampage and 
Finney ensembles, the coordination numbers were based 
on the position data. In the Rampage ensembles grains 
in a distance less than one grain diameter apart were 
considered touching, as the model allows small overlaps 
between grains. In the analysis of the Finney ensem- 
ble, grains were considered touching if the center to cen- 
ter distance r was less than 1.02 grain diameters. The 
choice of this distance was based on the average con- 
tact number obtained when different r's were assumed 
for touching neighbors: If grains were assumed touch- 
ing only when their center to center distance was 1 d or 
less, (fc) = 0.9. (fc) rapidly increased as the assumed dis- 
tance for touching grains increased, and at 1.02 d, the 
average coordination number was 6.72. Also, at 1.02 d, 
all grains (except for one) had at least three neighboring 
grains at distances closer or equal to r/d = 1.02, which is 
the stability criteria in 3D for grains shielded by granular 
bridges. 
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FIG. 10: The radial distribution functions for an ensemble 
of ductile grains (•, bold line) and for the Finney ensemble 
of dry, hard grains (thin line). The bold line was obtained 
by splines to the data. The inset shows the normalization 
function F{r), see Eq. II12II . (curve A) used to correct for 
the 3D finite size effects for the ductile system, in the region 
r/d e [0,5]. 




FIG. 11: The distribution of coordination numbers for the 
different ductile ensembles, as function of packing fraction. 
The curves A-F are splines to the data points of the packings 
at packing fractions c =0.64, 0.64, 0.66, 0.68, 0.70, and 0.75, 
respectively. 



E. Results (3D) 

To calculate RDFs, the center positions of the grains 
must be known, whereas in the 3D ductile experiments 
the positions of the top points of the grains were mea- 
sured. Although the top points are not the shifted posi- 
tions of the grain centers, they were used in calculating 
the RDF, as they represent correlations between specific 
points in the ensemble in the same way that the center 
positions do. Also, it should be directly comparable to 
the RDFs of hard granular ensembles, as in these the 
top positions are truly linear translations of the center 
positions. The detailed structure of an ensemble can 
only appear in the RDF if there is a sufficient amount 
of grains in the ensemble. Only the large ductile system 
contains enough grains to capture the main features of 
the RDF, thus unfortunately, we can only compare the 
ductile system at a packing fraction of c = 0.70 to the 
Finney ensemble (c — 0.636), and only at a resolution of 
dr — 0.1 d. The normalized RDF of the ductile ensemble 
is shown in Fig. ^| as is also the RDF of Finney's en- 
semble for dry granular media. From the figure, we see 
that the ductile ensemble is much less ordered than the 
hard granular ensemble, as all peaks and valleys in the 
RDF are small compared to the Finney RDF. We also 
note that the RDF of the ductile packing has a value for 
distances smaller than one grain diameter, which is not 
surprising, as this reflects the grains' ability to deform 
and thus obtain distances closer than one grain diameter. 
We see fewer clear peaks for the ductile ensemble than 
for the Finney ensemble, as deformed grains broaden the 
peaks. The split second peak in the RDF of the Finney 



ensemble is not present in the ductile ensemble. 

The coordination number distributions for the ductile 
systems are presented in Fig. The distributions shift 
to higher coordination numbers as the density increases, 
as is to be expected. One interesting observation is that 
two of the distributions for the small ensembles contain 
grains of a coordination number as low as three. This is 
possible for grains that are shielded by granular bridges, 
which are also found in the Finney and Rampage pack- 
ings to a similar degree (1.6% of the internal grains). Due 
to the cohesion between the ductile grains, the effect of 
granular bridging was expected to be more dominant in 
the looser ductile structures than in the Finney and Ram- 
page structures, which is not the case. For the denser 
ductile structures (c > 0.64) all grains have coordination 
numbers larger than three, but the width of the distri- 
butions remain the same within 4%. Figure El shows 
the mean coordination numbers as a function of packing 
fraction for the ensembles presented in this paper, and 
also for a few cases from the literature. The compaction 
of the ductile grain ensembles (filled circles) start in the 
neighborhood of the random close packed Finney ensem- 
ble (diamond), and then (fc) increases with c. The first 
five points on this curve were obtained from gravitation- 
ally compacted ensembles, thus their compaction process 
is the same. The last point, at c = 0.75, was obtained by 
mechanical compaction of a piston, and therefore does 
not necessarily represent the natural evolution of the 
gravitational ensembles. The numerical model (Ram- 
page) evolved along the curve represented by open cir- 
cles in Fig. 1121 during one simulation. It ends up also in 
the neighborhood of the Finney ensemble. Note that the 
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FIG. 12: Mean coordination numbers (fc) as a function of 
packing fraction c for several kinds of packings; The ductile 
ensembles (•), the Finney ensemble (o), Rampage ensembles 
(o), compacting bronze powder [l3l(*), compacted plasticine 
[2^ (a), hard ellipsoid ensemble 0] (+). The data for some 
crystalline structures of hard spheres are shown for reference: 
Simple cubic (□), Body centered cubic (□), and Face centered 
cubic (K). 



trajectory of the Rampage packing does not seem to be 
continued by the ductile ensembles, as the curves have 
different slopes at c = 0.64. Another ductile ensemble 
from the literature is presented in the figure (asterisks): 
Fischmeister et. al |l3j conducted compaction experi- 
ments on bronze powder by die compression at various 
pressures. The compacted ensemble was forced open af- 
ter compression, and the number of contacts and porosity 
measured. For the smallest packing fraction, the powder 
was sintered in order to mark the contact points of a 
non-compressed ensemble. The data at higher packing 
fractions were compacted at pressures of 0.2, 0.4, 0.6, 0.8 
and 1.0 GPa, respectively. The powder had an initial 
size distribution of 250-315 /im. The packing of hard el- 
lipsoids Q is marked as a plus sign in Fig. ^1 and the 
crystalline ensembles of simple cubic, body centered cu- 
bic and face centered cubic are plotted as various types 
of squares. 



IV. DISCUSSION 

That the shape of the grains is an important param- 
eter for the structure is evident from the experiment by 
Donev et al. , who studied the coordinations in an en- 
semble of ellipsoids and found a mean coordination num- 



ber of 9.82 in an ensemble of packing fraction c = 0.685. 
This value is much larger than that found in spherical 
ensembles at similar (but lower) c. The compaction of 
hard spherical grains to higher densities than the ran- 
dom close packed (rep) value (c = 0.63-0.64) |^ |2^ 
occurs by increasing the amount of crystalline regions. 
Pouliquen et al. 24] obtained a strongly crystalline pack- 
ing (c = 0.67) from mono-disperse glass spheres by hor- 
izontal shaking and slow filling of a container. However, 
handling hard spheres without fine-tuning the assembling 
procedure commonly results in values around c — 0.64. 
The compaction of deformable grains allows larger pack- 
ing fractions to be obtained than hard grain ensembles 
do, because the shape of individual grains can change. 
Whereas the crystallization of compacting hard grains 
would show up in the RDF, the evolution of the RDF 
during compaction is not given in ductile ensembles, as 
they might compact by structural ordering or grain de- 
formation. Shape changes from that of a sphere or a disk 
are expected to cause disordering (e.g., there are no crys- 
talline arrangement of pentagons). There are space fill- 
ing crystalline structures (squares, pyramids), but these 
are very unstable under any geometrical perturbation, as 
more than three edges meet at each vertex of the struc- 
ture. 

The radial distribution functions in two and three di- 
mensions as they were studied here, were both seen to 
lose some of the structure associated with hard grain ar- 
rangements as high packing fractions were approached. 
The RDF in 2D numerical hard disk ensembles has been 
shown [23 to develop its structure considerably through 
packing fractions c = 0.3-0.83. The split second, the 
third, and the fourth peaks were seen to develop be- 
tween c — 0.7-0.83, and the RDF at c = 0.83 resem- 
bles the RDF of the 2D ductile experiment in the ini- 
tial stage (Fig. EJ. Our 2D ductile system, when fur- 
ther compacted, gives less distinct splitting of peaks, and 
the RDF thus loses some of its detailed structure toward 
c = 0.99 (see Fig. E)). Interestingly, we saw that the peak 
at r/d = 1.695 corresponded to a cluster of four grains 
in the configuration shown in Fig. (Hjin the initial part of 
the experiment, the same configuration that contributes 
to the peak at r/d = -\/3 in hard disk ensembles. Thus, 
the initial part of our experiment is dominated by similar 
configurations as in hard disk ensembles, but with short- 
ened distances. At the end of our experiment, though, 
several different configurations contribute to this peak in 
the RDF, thus there are important structural changes 
compared to hard disk packings. 

A similar evolution is expected for 3D ensembles; A 
certain packing fraction must be reached for locally or- 
dered grain arrangements to occur, and for details in the 
RDF structure to emerge. When grain deformation be- 
comes important for further compaction, the RDF struc- 
ture becomes less distinct (Fig. [Tn|) . This 'destructuring 
should become influential in the 3D ductile grain ensem- 
bles at packing fractions close to the random close pack- 
ing, c ~ 0.64, assuming that frictional forces between the 
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grains prevent rearrangements. The RDF of the ductile 
3D system at c = 0.70 has no clear peaks after r/d — 3, 
and has much smaller and broader peaks than that of the 
2D system at c = 0.99. This difference could be due to 
the extra degree of freedom introduced by the third di- 
mension, or possibly the differences in rheology between 
the materials used in 2D and 3D, and how they react to 
the way they were compacted. Also a size distribution of 
the grains would tend to broaden the peaks in the RDF, 
but the size distribution is similar for our 2D and 3D 
experiments (3% and 2.3%, respectively). 

The evolution of the coordination number distribution 
suggests geometrical ordering to take place in the 2D duc- 
tile system. The distribution narrows considerably, and 
only the number of grains with six contacts increases, see 
Fig. |S1 Grains with only two contacts occur for a limited 
range of packing fractions c = 0.906-0.918, and these are 
due to shear motion, during which two neighboring or- 
dered regions align. In three dimensions, the evolution of 
the coordination number distribution is not as distinct. 
The average, (k), increases, but the width of the dis- 
tribution stays practically constant during compaction. 
The reason for this could be that the 3D ductile systems 
are not as compacted as the 2D system is, thus a marked 
narrowing of the coordination number distribution might 
not take place until larger packing fractions are reached. 
Also, the mean coordination number in 2D has a maxi- 
mum of six, while no established maximum exists in 3D. 
Thus, the distribution of k in 2D must narrow as the sys- 
tem gets denser, whereas the 3D system does not have 
this strict constraint. 

In 2D, the evolution of the average coordination num- 
ber (k) increased with packing fraction c toward the the- 
oretical value (fc) = 6, which applies to 2D space fill- 
ing structures (Voronoi cells). At the end of the 2D 
experiment, 1% of pore space remained, and according 
to Fig. the increase in average coordination number 
must be steep for c in the range 0.99-1. The evolution 
of (k) follows the straight line between the square lat- 
tice ((A;) = 4) and the space filling structure ((fe) = 6) 
up to c = 0.97. Although an increase of (fc) with c was 
expected, the linear evolution was not. The dominant 
crystalline structure is triangular from the start of the 
experiment, and regions of dense triangular structures 
(fc = 6) were thought to develop rapidly during the ini- 
tial compaction. 

3D deformable grains will approach a space filling 
structure of polyhedra with increasing packing fraction. 
This structure is not necessarily the packing's Voronoi 
structure, as the interfaces between neighboring grains 
then would have to be normal to the center-to-center 
distance. However, the statistical features of the Voronoi 
structure might be applicable at c = 1. No theoretical 
value exists for the average coordination at c = 1 in 3D, 
but a value around 13.6 was found by Bernal based on 
an experiment with compacted plasticine. The Voronoi 
tessellation of hard sphere packings at different packing 
fractions have shown that the average number of sides 



of a Voronoi cell, (/), decreases with increasing packing 
fraction [26L l27l | . From Voronoi analysis of the ductile 
ensemble at c = 0.70, we find that (/) = 14.5, and pre- 
sume that (/) in the ductile ensembles at c = 1 must be 
below this value. At a packing fraction of c = 1, all the 
faces of the polygon must be touching a neighbor, thus 
(/) = (fc) at c = 1. This suggests that had the ductile 
ensembles been allowed to compact to c = 1, their av- 
erage coordination number would be less than 14.5, and 
possibly close to Bernal's experimental value of 13.6. 

How the Voronoi cell structure changes during com- 
paction must depend on the compaction procedure. If 
the packing is compacted isomorphically, the relative po- 
sitions of the grains are fixed with respect to each other, 
and the Voronoi structure is constant. For grains on the 
square lattice, the Voronoi cells would remain squares 
throughout an isomorphic compaction, thus (fc) = 4 for 
all c. A random initial packing, or a non-isomorphic com- 
paction procedure, would instead result in an increase of 
the average coordination number toward (/). Thus the 
trajectories a compacting system makes in the ((fc), c)- 
plane is not fixed, but is expected to depend on initial 
configuration and compaction procedure. As seen from 
Fig. 1121 the ductile Play-Doh ensemble compacts in a 
self-compacting trajectory (apart from the last point on 
the curve), which evolves differently than the compacted 
bronze powder. The bronze powder was compacted by 
increasing uniaxial pressure, while our system sustained 
constant uniaxial pressure of gravity. Also, the Ram- 
page model of compaction has a trajectory leading into 
the region where the ductile compaction starts, but has 
a different slope at c = 0.64 than the trajectory of the 
ductile packings. The different compaction procedures 
might be an explanation for these different trajectories 
in the ((fc), c)-plane. 

Hard spherical ensembles can rather easily be com- 
pacted to the rep value, c ~ 0.64, and are well repre- 
sented by Finney's packing in the ((fc), c)-plane. This 
point might be fairly common for spherical ensembles, 
considering that both the Rampage model, and two of 
the ductile ensembles plot in the immediate neighbor- 
hood. Looser initial configurations might be obtained in 
packings of high cohesion or friction, which would change 
the starting point of the compaction in the {{k), c)-plane. 
That the bronze powder of Ref. [l^ plots at a higher (k) 
than the rep is thought to be due to the wide size distri- 
bution of grains, and the sintering procedure used. The 
hard elliptical grain ensemble [l^ was assembled in the 
same way as the Finney packing, but produced a much 
higher packing fraction and average coordination num- 
ber. The ability of the grains in the ductile packing to 
change their shape does not, however, take the ductile 
packing into the ellipsoid value of (fc) = 9.8 until c > 0.75 
for the Play-Doh ensembles, and c ~ 0.9 for the bronze 
powder (see Fig. 112(1 . Hence, the initial grain shape is 
an important factor for the initial compaction state in 
((fc), c)-plane. 
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V. CONCLUSION 

We have performed compactfon experiments on two- 
and three-dimensional packings of mono-dispcrsc ductile 
grains, and studied the evolution of their structure with 
increasing packing fraction. The radial distribution func- 
tion and coordination number distribution were found for 
both systems. The radial distribution function developed 
broader and smoother peaks, seemingly loosing ordered 
structure in both 2D and 3D. The 2D local configuration 
around grains in distances r/d = 1.695 ± 0.05 revealed 
that while mainly one configuration contributed to the 
corresponding peak in the initial RDF, two configura- 
tions contributed to the same peak at the most com- 
pacted stage of the experiment. The coordination num- 
ber distribution narrows considerably in the 2D-packing 
during the compaction, while no such narrowing is seen 
in the 3D-system. The average coordination number in 
3D evolves beyond that of hard spherical ensembles, and 
its evolution is discussed in relation to other granular en- 
sembles, as the compaction procedure and initial grain 



shape seems to be important parameters for the com- 
paction trajectory in the ((fc), c)-plane. 

This is a new approach to the study of compaction in 
granular ensembles, which emphasizes the importance of 
grain ductility on the evolving structure. Further work 
is necessary to understand the importance of compaction 
procedure, initial grain shape, and rheology for the final 
structure. In particular, 3D in-situ experimental inves- 
tigations of evolving geometry during compaction would 
give a strong contribution to the understanding of how 
local structure depends on compaction procedure. 
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